###############################################################################################################################################################################
#1. INTRODUCTION
#title: Analysis code (plotting of data) for Supran, Rahmstorf, and Oreskes (2022) 'Assessing ExxonMobil’s global warming projections'
#author: Geoffrey Supran
#date: 3 June 2022
#output(s): .csv in 'Temp Change/'; .pdf in 'Figures/'; .pdf in 'Figures/Graph Overlays/'
#description: This is the code used to generate the figures reported in Supran, Rahmstorf, and Oreskes (2022): 

# - Figure 1 (all panels except 1b): Historical temperature data on graphs scaled to original ExxonMobil graphs, enabling graphical overlays in figure 1 [section 5]
# - Figure 1 (panel 1b): Historical temperature data on 150,000-year timescale, scaled to original Black (1977) Exxon graph, enabling graphical overlay in figure 1 (panel 1b) [section 6]
# - Figure 2: Single time series ('spaghetti plot') of all projection data compared against averaged historical data [section 10]
# - Figure 3A: Temperature change per decade for projected versus historical data [section 3; average temperature changes also calculated and saved to dir = 'Results/Temp Change/']
# - Figure 3B: Implied transient climate response (iTCR) for projected versus historical data [section 4]
# - Figure S1A: Difference (projected minus observed) in temperature change per decade [section 8]
# - Figure S1B: Difference (projected minus observed) in iTCR [section 9]
# - Figure A of Print Summary: Historical CO2 data on graph scaled to Glaser (1982, fig.3) and Shaw (1984), enabling graphical overlay in figure A of Print Summary [section 7]

#These plots are constructed based on the data generated in Calculate.R and saved to dir = 'Data/' and dir = 'Results/'. Plots are output to dir = 'Results/Figures/'.
###############################################################################################################################################################################
#2. SETUP
library(readxl)
library(ggplot2)
library(egg)
library(tidyr)
library(forecast)
library(ggrepel)
library(boot)
setwd("###/")
###############################################################################################################################################################################
#3. MAIN CODE: Plot temperature change per decade for projected versus historical data
dir = 'Data/'
fname = 'TC_all_proj.rda'
load(file = paste(dir,fname,sep=''))
TC_all_proj2 = TC_all_proj

dir = 'Data/'
fname = 'TC_all_hist.rda'
load(file = paste(dir,fname,sep=''))
TC_all_hist2 = TC_all_hist

#Keep only rows of TC_all_hist2 that match TC_all_proj2 in terms of both category and CO2 scen
TC_all_hist2 = merge(TC_all_proj2, TC_all_hist2, by = c("cat", "CO2scen"))
TC_all_hist2 = TC_all_hist2[,-c(3:5)]
colnames(TC_all_hist2) = c('cat','CO2scen','TC','TC_high','TC_low')

TC_all_proj2$cat2 = paste('model',TC_all_proj2$CO2scen,sep = ' | ')
TC_all_hist2$cat2 = 'observation'
TC_all = rbind(TC_all_proj2, TC_all_hist2)
TC_all[,c(3:5)] = sapply(TC_all[,c(3:5)], as.numeric)
TC_all[,1] = as.factor(TC_all[,1])
TC_all[,2] = as.factor(TC_all[,2])
TC_all[,6] = as.factor(TC_all[,6])

##Calculate and include average TC of Exxon's projections and also of mainstream literature (from Hausfather et al. (2020)), for quotation in paper and plotting in figure 3A
#First calculate and include average/standard error of EXXON projections
f_med = function(data, indices){
  dt = data[indices,]
  c(mean(as.numeric(dt[,3])))
}

set.seed(1234)
bootstrap = boot(TC_all_proj2, f_med, R=10000)
TC_proj_avg = bootstrap$t0
TC_proj_se = sd(bootstrap$t)
temp = data.frame(cat = character(), CO2scen = character(), TC = numeric(), TC_high = numeric(), TC_low = numeric(), cat2 = character())
temp[1,1] = 'Averages'
temp[1,2] = 'nominal' #For plotting purposes, use the same factor nomenclature as for plotting other data of TC_all
temp[1,3] = TC_proj_avg
temp[1,4] = TC_proj_avg + 1.96*TC_proj_se
temp[1,5] = TC_proj_avg - 1.96*TC_proj_se
temp[1,6] = 'model | nominal' #For plotting purposes, use the same factor nomenclature as for plotting other data of TC_all
temp[,1] = as.factor(temp[,1])
temp[,2] = as.factor(temp[,2])
temp[,6] = as.factor(temp[,6])
TC_all = rbind(TC_all, temp)

#Second calculate and include average/standard error of mainstream literature projections (from Hausfather et al. (2020))
#Load results data from Hausfather et al. (2020)
dir = 'Raw data/'
df = read_excel(paste(dir,'2019 - Hausfather - model obs trends from SI.xlsx',sep=''), sheet = "main figures")
df2 = df[-c(69:72),c(1,3:6)] #Remove IPCC AR4 data not included in main analysis of Hausfather et al. (2020)
df3 = data.frame(df2[df2$variable == 'model_time',]) #Isolate 'model_time' data corresponding to temperature-versus-time metric of model projections

f_med = function(data, indices){
  dt = data[indices,]
  c(mean(dt[,3]))
}

set.seed(1234)
bootstrap = boot(df3, f_med, R=10000)
TC_Hausfather_avg = 10*bootstrap$t0 #Multiply by 10 to convert Hausfather et al. (2020)'s per-year data into per-decade 
TC_Hausfather_se = 10*sd(bootstrap$t) #Multiply by 10 to convert Hausfather et al. (2020)'s per-year data into per-decade

temp = data.frame(cat = character(), CO2scen = character(), TC = numeric(), TC_high = numeric(), TC_low = numeric(), cat2 = character())
temp[1,1] = 'Averages'
temp[1,2] = 'nominal' #For plotting purposes, use the same factor nomenclature as for plotting other data of TC_all
temp[1,3] = TC_Hausfather_avg
temp[1,4] = TC_Hausfather_avg + 1.96*TC_Hausfather_se
temp[1,5] = TC_Hausfather_avg - 1.96*TC_Hausfather_se
temp[1,6] = 'observation' #For plotting purposes, use the same factor nomenclature as for plotting other data of TC_all
temp[,1] = as.factor(temp[,1])
temp[,2] = as.factor(temp[,2])
temp[,6] = as.factor(temp[,6])
TC_all = rbind(TC_all, temp)

#Print summary spreadsheet of all TC values including averages. These averages are quoted in the main text: "We estimate the average warming projected by the 18 academic and government models to be 0.19 ± 0.03 °C per decade, which, within error, is the same as ExxonMobil’s average of 0.20 ± 0.04 °C per decade."
dir = 'Results/Temp Change/'
fname = 'TC_all.csv'
write.csv(TC_all, file = paste(dir,fname,sep=''))

#As a supplementary calculation, calculate average TC of Exxon's projections less those two that overlap with the 18 academic/government projections analyzed by Hausfather et al. (2020), for quotation in Supplementary Materials (SM): "[T]he average warming projected by ExxonMobil is unchanged after excluding the two overlaps: 0.20 ± 0.04 °C per decade."
models_Hausfather = c('1982 Weinberg; 1984 Callegari',
                      '2001 Albritton')

TC_all_proj2_nonHausfather = TC_all_proj2[!(TC_all_proj2$cat %in% models_Hausfather),]

f_med = function(data, indices){
  dt = data[indices,]
  c(mean(as.numeric(dt[,3])))
}

set.seed(1234)
bootstrap = boot(TC_all_proj2_nonHausfather, f_med, R=10000)
TC_proj_nonHausfather_avg = bootstrap$t0
TC_proj_nonHausfather_se = sd(bootstrap$t)
temp = data.frame(cat = character(), CO2scen = character(), TC = numeric(), TC_high = numeric(), TC_low = numeric(), cat2 = character())
temp[1,1] = 'Averages'
temp[1,2] = 'nominal' #For consistency with above tables, use the same factor nomenclature as for plotting other data of TC_all
temp[1,3] = TC_proj_nonHausfather_avg
temp[1,4] = TC_proj_nonHausfather_avg + 1.96*TC_proj_nonHausfather_se
temp[1,5] = TC_proj_nonHausfather_avg - 1.96*TC_proj_nonHausfather_se
temp[1,6] = 'model (less 2 projections overlapping with Hausfather et al.) | nominal'
temp[,1] = as.factor(temp[,1])
temp[,2] = as.factor(temp[,2])
temp[,6] = as.factor(temp[,6])
TC_all_SM = rbind(TC_all[c(nrow(TC_all)-1, nrow(TC_all)),], temp)
TC_all_SM$Bootstrapped_TwoSigma_StandardErrorOfMean = (TC_all_SM$TC_high - TC_all_SM$TC)/TC_all_SM$TC

#As another supplementary calculation, calculate average TC of EXXON'S OWN MODELS ONLY, for quotation in Table S1 of Supplementary Materials (SM): "Projections modeled by ExxonMobil scientists": "0.18 ± 0.04" °C per decade.
models_EM = c('1982 Glaser (Figure 3); 1984 Shaw',
              '1985 Flannery (Page 23)',
              '1985 Flannery (Page 24)',
              '1985 Hoffert (fig. 5.16)',
              '1994 Jain',
              '1997 Kheshgi',
              '2003 Kheshgi (fig. 7)',
              '2003 Kheshgi (fig. 8)')

TC_all_proj2_EM = TC_all_proj2[(TC_all_proj2$cat %in% models_EM),]

f_med = function(data, indices){
  dt = data[indices,]
  c(mean(as.numeric(dt[,3])))
}

set.seed(1234)
bootstrap = boot(TC_all_proj2_EM, f_med, R=10000)
TC_proj_EM_avg = bootstrap$t0
TC_proj_EM_se = sd(bootstrap$t)
temp = data.frame(cat = character(), CO2scen = character(), TC = numeric(), TC_high = numeric(), TC_low = numeric(), cat2 = character())
temp[1,1] = 'Averages'
temp[1,2] = 'nominal' #For consistency with above tables, use the same factor nomenclature as for plotting other data of TC_all
temp[1,3] = TC_proj_EM_avg
temp[1,4] = TC_proj_EM_avg + 1.96*TC_proj_EM_se
temp[1,5] = TC_proj_EM_avg - 1.96*TC_proj_EM_se
temp[1,6] = 'model (only models built/run by ExxonMobil) | nominal'
temp[,1] = as.factor(temp[,1])
temp[,2] = as.factor(temp[,2])
temp[,6] = as.factor(temp[,6])
TC_all_SM = rbind(TC_all_SM[,-ncol(TC_all_SM)], temp)
TC_all_SM$Bootstrapped_TwoSigma_StandardErrorOfMean = (TC_all_SM$TC_high - TC_all_SM$TC)/TC_all_SM$TC

#Print summary spreadsheet of all TC averages
dir = 'Results/Temp Change/'
fname = 'TC_averages_SM.csv'
write.csv(TC_all_SM, file = paste(dir,fname,sep=''))

#Plot figure 3A
cols=c(gray.colors(n=1, start = 0.5, end = 0.5), gray.colors(n=1, start = 0.5, end = 0.5), "black", "red")
dir = 'Results/Figures/'
plot = ggplot(data = TC_all, aes(x = cat2, y = TC, colour = cat2))+
  facet_grid(~cat,labeller = label_wrap_gen(width=20))+ #Width in labeller sets max number of characters before wrap
  geom_point() +
  geom_errorbar(aes(ymin = TC_low, ymax = TC_high), width=0.1, alpha = 1, linetype = "solid")+
  scale_color_manual(values=cols) +
  coord_cartesian(ylim = c(-0.01, 0.52)) + #Use coord_cartesian instead of scale_y_continuous in order that error bars display even if they exceed y-axis range.
  ylab("Temperature vs Time (degrees C/decade)") +
  theme_bw() %+replace%
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        legend.title = element_blank(),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=14),
        axis.title.y=element_text(size=14, angle = 90, margin = margin(t=0, r=15, b=0, l=0)),
        strip.text.x = element_text(size = 4, color = 'white'),
        strip.background = element_rect(colour = 'black', fill = "#1D71B8"),
        panel.grid.major.x = element_blank())
print(plot)
dev.off()

ggsave(filename = 'figure3A.pdf', plot = egg::set_panel_size(p=plot, width=unit(1.5, "cm"), height=unit(7, "cm")), device = "pdf", path = dir, scale = 1, width = 35, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#4. MAIN CODE: Plot iTCR for projected versus historical data
dir = 'Data/'
fname = 'TCR_all_proj.rda'
load(file = paste(dir,fname,sep=''))
TCR_all_proj2 = TCR_all_proj

dir = 'Data/'
fname = 'TCR_all_hist.rda'
load(file = paste(dir,fname,sep=''))
TCR_all_hist2 = TCR_all_hist

#Keep only rows of TCR_all_hist2 that match TCR_all_proj2 in terms of both category and CO2 scen
TCR_all_hist2 = merge(TCR_all_proj2, TCR_all_hist2, by = c("cat", "CO2scen"))
TCR_all_hist2 = TCR_all_hist2[,-c(3:5)]
colnames(TCR_all_hist2) = c('cat','CO2scen','TCR','TCR_high','TCR_low')

TCR_all_proj2$cat2 = paste('model',TCR_all_proj2$CO2scen,sep = ' | ')
TCR_all_hist2$cat2 = 'observation'
TCR_all = rbind(TCR_all_proj2, TCR_all_hist2)
TCR_all[,c(3:5)] = sapply(TCR_all[,c(3:5)], as.numeric)
TCR_all[,1] = as.factor(TCR_all[,1])
TCR_all[,2] = as.factor(TCR_all[,2])
TCR_all[,6] = as.factor(TCR_all[,6])

cols=c(gray.colors(n=1, start = 0.5, end = 0.5), gray.colors(n=1, start = 0.5, end = 0.5), "black", "red")
dir = 'Results/Figures/'
plot = ggplot(data = TCR_all, aes(x = cat2, y = TCR, colour = cat2))+
  facet_grid(~cat, labeller = label_wrap_gen(width=20))+ #Width in labeller sets max number of characters before wrap
  geom_point() +
  geom_errorbar(aes(ymin = TCR_low, ymax = TCR_high), width=0.1, alpha = 1, linetype = "solid")+
  scale_color_manual(values=cols) +
  coord_cartesian(ylim = c(-0.1, 5)) + #Use coord_cartesian instead of scale_y_continuous in order that error bars display even if they exceed y-axis range.
  ylab("Implied TCR (degrees C per 2x CO2)") +
  theme_bw() %+replace%
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        #panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        legend.title = element_blank(),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=14),
        axis.title.y=element_text(size=14, angle = 90, margin = margin(t=0, r=15, b=0, l=0)),
        strip.text.x = element_text(size = 4, color = 'white'),
        strip.background = element_rect(colour = 'black', fill = "#1D71B8"),
        panel.grid.major.x = element_blank())
print(plot)
dev.off()

ggsave(filename = 'figure3B.pdf', plot = egg::set_panel_size(p=plot, width=unit(1.5, "cm"), height=unit(7, "cm")), device = "pdf", path = dir, scale = 1, width = 35, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#5. MAIN CODE: Plot historical temperature data on graphs scaled to original ExxonMobil graphs, enabling graphical overlays in figure 1 (all panels except 1b)
dir = 'Data/'
fname = 'projections.rda'
load(file = paste(dir,fname,sep=''))
dim = startyears
dim$x0 = c(1850,1850,1960,1950,1850,1850,1840,1990,1800,1990,1990,1990)
dim$x1 = c(2050,2100,2100,2020,2100,2100,2100,2100,2200,2100,2100,2100)
dim$y0 = c(0,0,0,0,0,0,0,0,-5,0,0,0)
dim$y1 = c(2,2,1.2,0.6,2,2,2,2,10,2,2,2)
dim = dim[,-2]
dim = cbind(dim,offsets[,2:ncol(offsets)])

dir = 'Data/'
fname = 'historical.rda'
load(file = paste(dir,fname,sep=''))
hist$year = as.numeric(hist$year)

hist_avg = spread(hist, cat, dT)
temp = data.frame(hist_avg[,4:8])
hist_avg$'Historical (avg)' = rowMeans(temp, na.rm = T)
hist_avg2 = gather(hist_avg, cat, dT, -year, -sd_high, -sd_low)
hist = hist_avg2[,c(1,5,2:4)]

dir = 'Results/Figures/Graph Overlays/'
for(i in 1:nrow(dim)){
  
  #2003 Kheshgi is an anomalous figure in that the time period around which it is zeroed (1856-75) lies outside of the plot timespan (1990+). In this case, we therefore pre-offset the dataset so that it is zeroed around 1856-75, and then skip the zeroing procedure below (which only works within the timespan of the plotted region).
  if(dim$cat[i] == '2003 Kheshgi (fig. 7)' | dim$cat[i] == '2003 Kheshgi (fig. 8)'){
    hist_temp = hist[hist$cat == 'Historical (avg)',]
    hist_temp2 = data.frame(dT = hist_temp[hist_temp$year>=1856 & hist_temp$year<=1875,2]) #Isolate portion of historical temperature dataseries to be averaged for comparison against preindustrial range of projection dataseries
    offset_hist = mean(hist_temp2$dT)
    hist$dT = hist$dT - offset_hist
  }
  
  hist2 = hist[hist$year >= dim[i,2] & hist$year <= dim[i,3],] #Select data only over x-axis range of relevant plot
  
  #For averaged historical dataseries, calculate its average value over the same preindustrial period as projection data, then add the difference between that average and the average of projection's preindustrial series to the historical series in order that their averages are equal over the same preindustrial range.
  hist_temp = hist2[hist2$cat == 'Historical (avg)',]
  hist_temp2 = lowess(hist_temp$year, hist_temp$dT, f = 0.3) #LOWESS smoothing with span of f = 0.3
  hist_temp2 = data.frame(year = hist_temp2$x, dT = hist_temp2$y)
  
  if(dim$cat[i] != '2003 Kheshgi (fig. 7)' && dim$cat[i] != '2003 Kheshgi (fig. 8)'){ #In the anomalous case of 2003 Kheshgi, skip this zeroing step (it has been zeroed above).
    hist_temp3 = data.frame(dT = hist_temp2[hist_temp2$year>=dim[i,6] & hist_temp2$year<=dim[i,7],2]) #Isolate portion of historical temperature dataseries to be averaged for comparison against preindustrial range of projection dataseries
    offset_hist = mean(hist_temp3$dT)
    
    diff = offsets[i,4] - offset_hist
    hist_temp2$dT = hist_temp2$dT + diff
  }
  hist3 = hist_temp2
  
  #Crop the data range of select series in order that overlaid plotted data fits within axes of original EM plots
  if(dim$cat[i] == '1982 Glaser (Figure 3); 1984 Shaw'){
    hist3 = hist3[hist3$year >= 1979,]
  }
  if(dim$cat[i] == '1982 Weinberg; 1984 Callegari'){
    hist3 = hist3[hist3$dT<=0.8,]
  }
  
  plot = ggplot(data = hist3, aes(x = year, y = dT))+
    geom_line(colour = "red")+
    coord_fixed(0.8) +
    scale_x_continuous(limits = c(1800, 2100), breaks = c(dim[i,2], dim[i,3])) +
    scale_y_continuous(limits = c(-2,2), breaks = c(dim[i,4], dim[i,5]))+
    theme_bw() +
    theme(panel.grid.major = element_blank(), 
          panel.grid.minor = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA, size = 1),
          axis.ticks.length = unit(-3, "mm"),
          axis.text.y.right = element_blank(),
          axis.text.x.top = element_blank(),
          axis.title.x.bottom = element_blank(),
          axis.title.y.right = element_blank(),
          axis.title.x.top = element_blank(),
          legend.title = element_blank(),
          axis.text.x = element_text(margin = margin(t = 0.5, unit = "cm"), size=16),
          axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=16),
          axis.title.y=element_text(size=16, margin = margin(t=0, r=15, b=0, l=0))) + 
    ylab("Temperature Change (degrees C)")
  print(plot)
  dev.off()
  
  ggsave(filename = paste('OVERLAY_',toString(dim[i,1]),'.pdf',sep=''), plot = egg::set_panel_size(p=plot, width=unit(14, "cm"), height=unit(11, "cm")), device = "pdf", path = dir, scale = 1, width = 30, height = 25, units = c("cm"),
         dpi = 300)
}
###############################################################################################################################################################################
#6. MAIN CODE: Plot historical temperature data on 150,000-year timescale, scaled to original Black (1977) Exxon graph, enabling graphical overlay in figure 1 (panel 1b)
dir = 'Raw data/'
df = read_excel(paste(dir,'glacial.xlsx',sep=''), sheet = "Sheet1")

df_long = data.frame(df[,c(1,10)])
colnames(df_long) = c('year','dT')

df_long2 = data.frame(lowess(df_long$year, df_long$dT, f = 0.001)) #LOWESS smoothing with span of f = 0.3
colnames(df_long2) = c('year','dT')

dir = 'Results/Figures/Graph Overlays/'
plot = ggplot(data = df_long2, aes(x = year, y = dT))+
  geom_line(colour = "red")+
  coord_fixed(0.8) +
  scale_x_continuous(limits = c(-150, 25), breaks = c(-150,25)) +
  scale_y_continuous(limits = c(-6,3), breaks = c(-6,-4,-2,-2.7777777,0,2,2.7777777))+
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        axis.text.y.right = element_blank(),
        axis.text.x.top = element_blank(),
        axis.title.x.bottom = element_blank(),
        axis.title.y.right = element_blank(),
        axis.title.x.top = element_blank(),
        legend.title = element_blank(),
        axis.text.x = element_text(margin = margin(t = 0.5, unit = "cm"), size=16),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=16),
        axis.title.y=element_text(size=16, margin = margin(t=0, r=15, b=0, l=0))) + 
  ylab("Temperature Change (degrees C)")
print(plot)
dev.off()

ggsave(filename = paste('OVERLAY_1977 Black_LongTimeScale.pdf',sep=''), plot = egg::set_panel_size(p=plot, width=unit(14, "cm"), height=unit(11, "cm")), device = "pdf", path = dir, scale = 1, width = 30, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#7. MAIN CODE: Plot historical CO2 data on graph scaled to Glaser (1982, fig.3) and Shaw (1984), enabling graphical overlay in figure A of Print Summary
dir = 'Raw data/'
df = read_excel(paste(dir,'co2_annmean_mlo.xlsx',sep=''), sheet = "co2_annmean_mlo")
df2 = df[c(56:nrow(df)),c(1,2)]
colnames(df2) = c('year','CO2')
df2[,c(1:2)] = sapply(df2[,c(1:2)], as.numeric)

dir = 'Results/Figures/Graph Overlays/'
plot = ggplot(data = df2, aes(x = year, y = CO2))+
  geom_line(colour = "red")+
  coord_fixed(0.8) +
  scale_x_continuous(limits = c(1960, 2100), breaks = c(1960,2100)) +
  scale_y_continuous(limits = c(300,660), breaks = c(300,580))+
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        axis.text.y.right = element_blank(),
        axis.text.x.top = element_blank(),
        axis.title.x.bottom = element_blank(),
        axis.title.y.right = element_blank(),
        axis.title.x.top = element_blank(),
        legend.title = element_blank(),
        axis.text.x = element_text(margin = margin(t = 0.5, unit = "cm"), size=16),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=16),
        axis.title.y=element_text(size=16, margin = margin(t=0, r=15, b=0, l=0))) + 
  ylab("CO2 (ppm)")
print(plot)
dev.off()

ggsave(filename = paste('OVERLAY_1982 Glaser (Figure 3); 1984 Shaw_CO2.pdf',sep=''), plot = egg::set_panel_size(p=plot, width=unit(14, "cm"), height=unit(11, "cm")), device = "pdf", path = dir, scale = 1, width = 30, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#8. MAIN CODE: Plot difference (projected minus observed) in temperature change per decade
dir = 'Data/'
fname = 'TC_all_diff.rda'
load(file = paste(dir,fname,sep=''))

TC_all_diff$cat2 = paste('model',TC_all_diff$CO2scen,sep = ' | ')
TC_all_diff[,c(3:5)] = sapply(TC_all_diff[,c(3:5)], as.numeric)
TC_all_diff[,1] = as.factor(TC_all_diff[,1])
TC_all_diff[,2] = as.factor(TC_all_diff[,2])
TC_all_diff[,6] = as.factor(TC_all_diff[,6])

cols=c(gray.colors(n=1, start = 0.5, end = 0.5), gray.colors(n=1, start = 0.5, end = 0.5), "black")
dir = 'Results/Figures/'
plot = ggplot(data = TC_all_diff, aes(x = cat2, y = TCdiff, colour = cat2))+
  #facet_grid(~cat, labeller = label_wrap_gen(width=20))+ #Width in labeller sets max number of characters before wrap
  facet_grid(~cat,labeller = label_wrap_gen(width=20))+ #Width in labeller sets max number of characters before wrap
  geom_point() +
  geom_errorbar(aes(ymin = TCdiff_low, ymax = TCdiff_high), width=0.1, alpha = 1, linetype = "solid")+
  scale_color_manual(values=cols) +
  coord_cartesian(ylim = c(-0.3, 0.3)) + #Use coord_cartesian instead of scale_y_continuous in order that error bars display even if they exceed y-axis range.
  ylab("Difference (projection vs observation) in Temperature vs Time (degrees C/decade)") +
  geom_hline(yintercept=0, color = "red", linetype = 'dashed', size=1) +
  theme_bw() %+replace%
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        legend.title = element_blank(),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=14),
        axis.title.y=element_text(size=14, angle = 90, margin = margin(t=0, r=15, b=0, l=0)),
        strip.text.x = element_text(size = 4, color = 'white'),
        strip.background = element_rect(colour = 'black', fill = "#1D71B8"),
        panel.grid.major.x = element_blank())
print(plot)
dev.off()

ggsave(filename = 'figureS1A.pdf', plot = egg::set_panel_size(p=plot, width=unit(1.5, "cm"), height=unit(7, "cm")), device = "pdf", path = dir, scale = 1, width = 35, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#9. MAIN CODE: Plot difference (projected minus observed) in iTCR
dir = 'Data/'
fname = 'TCR_all_diff.rda'
load(file = paste(dir,fname,sep=''))

TCR_all_diff$cat2 = paste('model',TCR_all_diff$CO2scen,sep = ' | ')
TCR_all_diff[,c(3:5)] = sapply(TCR_all_diff[,c(3:5)], as.numeric)
TCR_all_diff[,1] = as.factor(TCR_all_diff[,1])
TCR_all_diff[,2] = as.factor(TCR_all_diff[,2])
TCR_all_diff[,6] = as.factor(TCR_all_diff[,6])

cols=c(gray.colors(n=1, start = 0.5, end = 0.5), gray.colors(n=1, start = 0.5, end = 0.5), "black")
dir = 'Results/Figures/'
plot = ggplot(data = TCR_all_diff, aes(x = cat2, y = TCRdiff, colour = cat2))+
  facet_grid(~cat,labeller = label_wrap_gen(width=20))+ #Width in labeller sets max number of characters before wrap
  geom_point() +
  geom_errorbar(aes(ymin = TCRdiff_low, ymax = TCRdiff_high), width=0.1, alpha = 1, linetype = "solid")+
  scale_color_manual(values=cols) +
  coord_cartesian(ylim = c(-2, 3)) + #Use coord_cartesian instead of scale_y_continuous in order that error bars display even if they exceed y-axis range.
  ylab("Difference (projection vs observation) in Implied TCR (degrees C per 2x CO2)") +
  geom_hline(yintercept=0, color = "red", linetype = 'dashed', size=1) +
  theme_bw() %+replace%
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        legend.title = element_blank(),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=14),
        axis.title.y=element_text(size=14, angle = 90, margin = margin(t=0, r=15, b=0, l=0)),
        strip.text.x = element_text(size = 4, color = 'white'),
        strip.background = element_rect(colour = 'black', fill = "#1D71B8"),
        panel.grid.major.x = element_blank())
print(plot)
dev.off()

ggsave(filename = 'figureS1B.pdf', plot = egg::set_panel_size(p=plot, width=unit(1.5, "cm"), height=unit(7, "cm")), device = "pdf", path = dir, scale = 1, width = 35, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################
#10. MAIN CODE: Plot single time series ('spaghetti plot') of all projection data compared against averaged historical data
dir = 'Data/'
fname = 'projections.rda'
load(file = paste(dir,fname,sep=''))

dir = 'Data/'
fname = 'historical.rda'
load(file = paste(dir,fname,sep=''))

hist$year = as.numeric(hist$year)
hist_avg = spread(hist, cat, dT)
temp = data.frame(hist_avg[,4:8])
hist_avg$'Historical (avg)' = rowMeans(temp, na.rm = T)
hist_avg2 = gather(hist_avg, cat, dT, -year, -sd_high, -sd_low)
hist = hist_avg2[,c(1,5,2:4)]
hist2 = hist[hist$cat == 'Historical (avg)',]
hist2 = data.frame(hist2[,-c(3:4)])

hist3 = lowess(hist2$year, hist2$dT, f = 0.3) #LOWESS smoothing with span of f = 0.3
hist3 = data.frame(year = hist3$x, dT = hist3$y)
hist3$dT = hist3$dT - hist3[hist3$year == 1900,'dT'] #Zero historical data at first year (1900) plotted in Fig. 2.

hist3$cat = 'historical'
hist3$cat2 = 'Observations'
hist3$label = NA

all_dT_nominal = data.frame()
cats = unique(proj$cat)
for(i in 1:length(cats)){
  temp = proj[proj$cat == cats[i],]
  temp2 = temp[temp$CO2scen == 'nominal',]
  
  if(cats[i] == '1997 Kheshgi'){ #For the unique case of 1997 Kheshgi, which has only two data points, the earlier of which (1995) is two years before startyears[i,2] (1997), we retain the 1995 data point for plotting purposes.
    temp3 = temp2
    offset = temp3$dT[1] - hist3[hist3$year==temp2$year[1],2] #Tie each projection data set to historical data at the start year of the projection data set, by computing and then removing the offset between projected and historical data on the start year.
  }else{ 
    temp3 = temp2[temp2$year >= startyears[i,2],]
    offset = temp3$dT[1] - hist3[hist3$year==startyears[i,2],2] #Tie each projection data set to historical data at the start year of the projection data set, by computing and then removing the offset between projected and historical data on the start year.
  }  
  
  temp3$dT = temp3$dT - offset
  temp3$cat2 = 'Projections'
  temp3$label = NA
  val = temp3$cat[which(temp3$year == max(temp3$year))]
  ind = which(cats %in% val)
  temp3$label[which(temp3$year == max(temp3$year))] = ind
  
  all_dT_nominal = rbind(all_dT_nominal, temp3)
  #  }
}

all_dT_nominal2 = all_dT_nominal[,-c(3:6)]

df = rbind(hist3,all_dT_nominal2)
df$cat = factor(df$cat, levels=unique(df$cat))

cols=c("red", gray.colors(n=length(unique(all_dT_nominal2$cat)), start = 0.9, end = 0.1))
sizes = c(3,0.5)
dir = 'Results/Figures/test/'
plot = ggplot(data = df, aes(x = year, y = dT, col = cat))+
  geom_line(aes(size = cat2))+
  scale_size_manual(values = sizes) +
  coord_fixed(0.8) +
  scale_x_continuous(limits = c(1900, 2100), sec.axis = dup_axis(), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.2,3.5), sec.axis = dup_axis(), expand = c(0, 0))+
  theme_bw() +
  scale_color_manual(values=cols) +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, size = 1),
        axis.ticks.length = unit(-3, "mm"),
        axis.text.y.right = element_blank(),
        axis.text.x.top = element_blank(),
        axis.title.x.bottom = element_blank(),
        axis.title.y.right = element_blank(),
        axis.title.x.top = element_blank(),
        legend.title = element_blank(),
        axis.text.x = element_text(margin = margin(t = 0.5, unit = "cm"), size=16),
        axis.text.y = element_text(margin = margin(r = 0.5, unit = "cm"), size=16),
        axis.title.y=element_text(size=16, margin = margin(t=0, r=15, b=0, l=0))) + 
  ylab("Temperature Change (degrees C)")
print(plot)
dev.off()

ggsave(filename = 'figure2.pdf', plot = egg::set_panel_size(p=plot, width=unit(15, "cm"), height=unit(12, "cm")), device = "pdf", path = dir, scale = 1, width = 35, height = 25, units = c("cm"),
       dpi = 300)
###############################################################################################################################################################################